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ABSTRACT 

We present models of turbulent mixing at the boundaries between hot (T ~ 10 6-7 K) and warm 
material (T ~ 10 4 K) in the interstellar medium, using a three-dimensional magnetohydrodynamical 
code, with radiative cooling. The source of turbulence in our simulations is a Kelvin-Hclmholtz 
instability, produced by shear between the two media. We found, that because the growth rate of the 
large scale modes in the instability is rather slow, it takes a significant amount of time (~ 1 Myr) 
for turbulence to produce effective mixing. We find that the total column densities of the highly 
ionized species (C IV, N V, and O VI) per interface (assuming ionization equilibrium) are similar to 
previous steady-state non-equilibrium ionization models, but grow slowly from log N ~ 10 11 to a few 
x 10 12 cm -2 as the interface evolves. However, the column density ratios can differ significantly from 
previous estimates, with an order of magnitude variation in N(C IV)/N(0 VI) as the mixing develops. 
Subject headings: ISM: general — ISM: structure — magnetohydrodynamics: MHD — turbulence 



1. INTRODUCTION 

The interstellar medium (ISM) is a very complex en- 
tity. It is extremely rich in structure, highly turbu- 
lent, and embedded in a dynamically important mag- 
netic fi eld. Altho ugh the concept of "phases" is up for 
debate (|Cox l2005D , we can safely say that the ISM shows 
regions of distinct physical conditions, which range from 
cold and molecular, to hot and ionized. A factor that 
controls the structure of the ISM to a large extent, is the 
balance between heating and cooling processes. In this 
work, we are mainly concerned with the interfaces be- 
tween the hot (T ~ 10 6 K) and warm (T - 10 4 K) media, 
as well as the transfer of heat between them. At temper- 
atures in between the hot and war m media the efficiency 
of rad iative cooling is maximal (see lSutherland fc Dopital 
1993). Material at such temperatures cools very quickly, 
therefore should be rarely observed. However, absorption 
line studies i n the far ultraviolet (FUV) have found oth- 
erwise (e. sr. iSavasre et al.ll20 03^. suggesting substantial 
amounts of plasma at temperatures of T ~ 10 5 K. 

It was realized by Bcgchnan & Fabian ( 1990, hereafter 
BF90), that the interfaces between hot and warm media 
might well be dominated by "turbulent mixing layers" . 
The general idea was that turbulence produces an ex- 
change of material and energy at the boundary between 
the two media, providing a steady supply of plasma at 
T ~ 10 5 K, balancing the l osses d ue to radiation. Later, 
ISlavin. Shull fc BegelmaJ IpM hereafter SSB93) pre- 
dicted that line ratios of highly ionized species (C IV, 
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N V, O VI, and Si IV) differ significantly from those 
in purely photoionized gas, radiative shocks, and con- 
duction fronts, thus providing a useful diagnostic tool 
of the physical conditions in turbulent mixing layers. 
Their model however, relied on several simplifying as- 
sumptions. For instance they characterized the mixing 
layer with a single mean temperature T, regardless of 
the position in the mixing layer. The model is one- 
dimensional, unmagnetized, and parameters were chosen 
to correspond observations. The obvious requirement for 
this type of heat transport is the development of turbu- 
lence at the boundary of the two media. In the ISM there 
are many outflows and plasma instabilities that can lead 
to this situation. For example, in a supernova explo- 
sion Rayleigh- Taylor instabilities occur when the hot and 
low-density ejecta tries to accelerate a colder and denser 
medium. This produces a boundary layer that breaks 
up into small (colder) high density clumps left inside a 
hot cavity. Another physical process that can lead to a 
turbulent mixing layer is the Kelvin-Helmholtz (K-H) in- 
stability, which can occur when shear is present between 
two fluids. This is possible, for instance, at the edges of 
high velocity cl ouds (HVCs) moving through t he galactic 
hot corona fsee lWakker fc van Woerdenlll997lL 

However, the presence of a magnetic field can influence 
the dynamics of the K-H instability. It is well known 
that flows with an Alfven speed less th an the shear ve- 
locity can become stable (see for instance |Chandrasekharl 
119611) . A detail ed numerical study of the K-H instability 
can be found in lFrank et al.l |)1996|) : IRvu. Jones fc Frankl 
(2000). Enhancement of thermal difussivity in magne- 
tized plasmas du e to turbulent motions is discussed in 
iCho et alJ pOO^ . 

In this work, we study the formation and evolution of 
turbulent mixing layers by means of the K-H instability. 
We use a magnetohydrodynamic (MHD) code, which in- 
cludes radiative cooling. The layout of the paper is as 
follows: in §2 we provide with a brief review of the the- 
ory behind turbulent mixing layers, in §3 we describe the 
code and our numerical setup. The results, including es- 
timates of column densities and line ratios of highly ion- 
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ized species can be found in §4, followed by a summary 
in §5. 

2. TURBULENT MIXING LAYERS 

The idea of a turbulent mixing layer was proposed in 
BF90 as an important heat transfer mechanism between 
two media at different temperatures. The basic picture 
proposed is that turbulence at the boundary of the two 
fluids will provide a constant input to the intermediate 
temperature mixture. If energy is conserved, this mix- 
ing layer would grow indefinitely, and eventually all of 
the material will be at such intermediate temperature. 
To prevent this they proposed a steady state in which 
the energy lost by radiation (radiative cooling is most 
efficient precisely at such intermediate temperatures) is 
balanced by a turbulent heat flux into the mixing layer. 
As it is usual, most of the energy in the turbulence is 
on the largest scales, and cascaded down to a dissipative 
level. Thus the model of BF90 is basically a three-phase 
steady state fluid in which the losses due to radiation 
are balanced with energy entrained by turbulence into 
the intermediate temperature zone. In their model, the 
temperature of the mixing layer is determined by mass 
flux balance: 



T 



thhotThot + rricoidXcoid 
rhhot + rhcoid 



(1) 



where, rn-hot and m co \d are the mass flux rates from the 
hot and from the cold phase into the layer, That and T co id 
are the temperatures of the hot and cold phases respec- 
tively. Since neither the heat transfer nor the mixing will 
be perfect, BF90 introduced two efficiency factors: r\hot, 
the fraction of mass and energy that is deposited by the 
hot medium to the mixing layer, and r\ C oid, an efficiency 
for the hydrodynamic mixing. With these factors the 
mass flux rates become 



mhot = r]hotPhotVti 

rncold = Vcold (photPc 



,1/2 



(2) 
(3) 



and the resulting intermediate temperature will be: 



Tf 



Vhot + Vcold {Tcold/Thot} 1 ^ 2 
Vcold + Vhot (Tcold/Thot) 1 ^ 2 
£ (TcoldThot) 1 ^ ■ 



(T C old Thot) 



1/2 



(4) 

The definition of the two different efficiencies reveals 
what the authors had in mind as the mechanism that 
provides the mixing. The efficiency associated with the 
entrainment of hot gas simply corresponds to an enthalpy 
flux ^rjhpvt, where p is the pressure and vt the turbulent 
velocity. This yields eq.J2J), where ph is the mass den- 
sity of the hot medium and the contribution of turbu- 
lent kinetic energy has been neglected. The cold gas is 
then pulled into the hot medium by turbulent eddies at a 
rate rj c p c l c /t(l c ). Turbulent eddies are formed on scales 
l c < It, where It is the total thickness of the layer, and 
t(l c ) can be thought as an eddy turnover time for eddies 
of size l c . BF90 used the Kelvin- Hclmholtz growth rate 

tKH(lc) ~ (pc/ PhY^lc/vt, where p c is the mass density 
of the cold medium, leading to eq.©. However, the form 

of the timescale ~ (p c /p/s) W w < ^ s n °t exclusive for 
the K-H instability but suitable for fully developed tur- 
bulence in general. The main uncertainty in the BF90 



model lies in these efficiency factors. It is also derived 
assuming fully developed turbulence (i.e. rapid mixing), 
and therefore does not include the effects of cooling to 
the dynamical development of turbulence. 

Later, SSB93 expanded on the ideas of BF90 and ran a 
grid of models based on one-dimensional, instantaneous, 
steady-state mixing, characterized principally by T, and 
v t . These included the effects of non-equilibrium ion- 
ization and self-photoionization of the gas in the mixing 
layer, but adopted a somewhat ad hoc value for £ and 
the efficiencies (rj's). 

In this work we focus on the dynamical formation and 
development of the mixing layer. We do not include ef- 
fects of non-equilibrium or self-photoionization, but we 
measure the the actual T, which is result a continu- 
ous distribution of temperatures that range from Thot 
to T co id, e instead of fixing any efficiency factor. 

3. OUR MODEL 
3.1. The code 
We solve the following system of equations: 

dp 



at 



+ V • (pv) = 0, 



<9v _ 1 (VxB)xB „ 

— + V • Vv + -Vp - ^ -r 1 = 0, 

ot p 4irp 



dp 

at 



v • Vp + 7pV • v = 0, 



<9B 



- V x (v x B) = 0, 



(5) 
(6) 
(7) 
(8) 



with V • B = 0. Here p is the mass density, v is the 
velocity, p is the pressure, and B the magnetic field. 
We use an ideal gas equation of state p = (7 — l)pu, 
where 7 is the ratio of the specific heats (7 = c p /c v ), 
and u is the specific internal energy. We solve equations 
(Ell-dHt using a MUSCL-type scheme with HLL fluxes 
ijHarten. Lax fc jyjm_ijeialLL9_8_3j) and monotoniz ed central 
limiter fsee lKurganov. Noelle fc Petroval2 001'l. on a reg- 
ular Cartesian grid. We use the the flux-inter polated con- 
strained transport (or "flux-CT") scheme (|X6thl |200C() 
to enforce the divergence free condition of the mag- 
ne tic field. A similar version of the c ode was used 
bv lGammie . McKi nnev. J., fc Tothl l)2003j) for relativistic 
MHD calculations. Our code gives satisfactor y results for 
standard shock tube tests (see, for example iBrio fc Wul 
119881 IRvu fe Jones! 1 1995]) . When compared with our 
earlier hybrid ENO c ode (see an isothermal version in 
ICho fc Lazari a*nl2002jL our current code runs faster, and 
yet gives consistent results. The overall scheme is second- 
order accurate. After updating the system of equations 
along xl direction, we repeat similar procedures for x2 
and x3 directions, with the appropriate rotation of in- 
dexes. We use a two-stage Runge-Kutta method for 
time integration. The cooling at a given temperature 
is assumed to be of that of plasma in collisional ion- 
ization equilibrium, with solar metallicity, as given in 
iBeni amm. Benson fc Coxl (|2001^l . We consider 7 = 5/3, 
and solve the time evolution of specific internal energy. 
For the runs with cooling, this is applied before updating 

6 Instead of T co id we will call it T marm because our choice of 
parameters are to coincide with the warm phase of the ISM (at 
T ~ 10 4 K). 
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equations JSJ-JSJ, using an implicit scheme (i.e. we sim- 
ply subtract the internal energy that will be lost in the 
time-step St). We do not explicitly include heating due 
to absorption of radiation or thermal conduction. The 
time-step for the MHD part is determined by the stan- 
dard "Courant" condition. When we include cooling, we 
monitor the minimal cooling time in the computational 
box and compare it with that of the MHD part. If the 
minimum cooling time is shorter, it replaces the time- 
step. However, this was very rarely the case, and typ- 
ically the dynamical time was shorter than the cooling 
time (by one or two orders of magnitude in most cases). 

We do not include thermal conduction explicitly. The 
diffusion of heat at small scales is determined in our 
simulations by the numerical diffusion, whose proper- 
ties may be different from those of diffusion in the ac- 
tual interstellar gas. However, this does not greatly af- 
fect our final results if the heat transport is determined 
by turbulence. A discussion about turbulent and ther- 
m al conduct i on in a magnetized medium can be found 
in iLazarianl (J2006) . Turbulence provides an effective 
diffusion coefficient ~ l/3wt ur hLj nj - , where v tur b is the 
turbulent velocity, and Li n j is the scale of energy in- 
j ection, correspo nding to the size of the largest eddies 
(ICho et alJl2003h . We will see that in our models the 
injection scale is utterly determined by the size of our 
computational box. If the turbulent diffusion coefficient 
is much larger than the thermal diffusion coefficient the 
heat transfer is dominated by turbulence. An impor- 
tant property of turbulent heat transfer, as well as the 
other transport processes, is that they do not depend 
on the microscopic diffusivity. Indeed, provided that the 
gas is turbulent, the mixing and the heat transport are 
happening approximately over one large eddy turnover 
time. If the diffusivity on the atomic level decreases, 
the turbulent cascade goes to yet smaller scales ensur- 
ing heat transport that still scales as above. For typi- 
cal parameters of turbulent mixing layers in our models 



(that is T ~ 10 5 K, n 



10- 



cm 



10 pc, 



Vturb ~ 20 kms 1), one obtains a Spitzer thermal diffu- 
sion coefficient of ks p < 10 24 cm 2 s" 1 and a turbulent 
diffusion coefficient of Kturb > 10 25 cm 2 s _1 . The differ- 
ence is modest, however the presence of a magnetic field 
will further suppress thermal conduction. Moreover, the 
magnetic field does not have to be dynamically dominant 
to produce an important effect on the thermal conduc- 
tivity. And, in some sense, our unmagnetized models 
are equivalent to models with a very weak magnetic field 
(ubiquitous in the ISM), that is sufficient to suppress 
thermal conductivity significantly. At the beginning of 
the simulations turbulence is not fully developed and the 
turbulent diffusivity is small, but at the same time the 
magnetic field is aligned with the interface, dramatically 
decreasing electron conduction as well. At more evolved 
stages however, turbulence will develop and thermal con- 
duction will be less suppressed as the magnetic field be- 
comes entangled. Estimates by iNaravan fc Medvedevl 
l)2001[) suggest that for fully developed turbulence the 
thermal conductivity is decreased by a factor of ~ 5 from 
the Spitzer value. Their model is, however, rather re- 
strictive as only turbulence with Alfven Mach number 
(Ma) equal to unity is discussed. For both Ma much 
larger, and much smaller than unity the electron thermal 




Fig. 1. — Schematic of the calculation cube, with hot gas moving 
upward on the left, and initially static warm gas to the right. 



conductivity is less fsee lLazaria n1 l200l . 

3.2. The numerical setup 

We start with warm gas, at a temperature of T warm — 
10 4 K, with a relatively high hydrogen number density 
n wa rm = 0.1 cm -3 , and no mean motion, at the right side 
of the computational box. At the opposite side we have a 
low-density (rihot), hot medium (Thot), moving upward (z 
direction) with a velocity Vhot- In some of our simulations 
we include a magnetic field B z warm threading the warm 
gas aligned in the vertical direction, with a magnitude 
corresponding to a /3 warm ~ Pgas/Pmag ~ 1. A sketch of 
the computational domain is presented in Figure ^ 

The two media are initially in total pressure equilib- 
rium. The transition region (between hot and warm me- 
dia) follows a hyperbolic tangent profile in the x direc- 
tion, that initially occupies ~ l/10th. of the box size. 
The division line for the runs in a 144 3 grid was set at 
the middle of the box. The vertical direction in all cases 
correspond to a physical scale of L z = 10 pc. In the cases 
without cooling, the simulations can be rescaled arbitrar- 
ily, but when cooling is included the scale lengths are 
fixed by the value of the cooling rate used. The boundary 
conditions are periodic in the z and y directions. For the 
x direction in the warm side of the box (right) we have 
a reflective boundary, and on the hot side (left) we have 
a source boundary condition. The latter reinforces the 
initial condition after each time-step, acting as a reser- 
voir of hot material, that helps to balance the energy 
lost trough radiation. We initialized the computational 
cube with sinusoidal perturbations (32 harmonics, with 
random phases) in the component of the velocity normal 
to the interface of the two media (v x ). The shear at the 
boundary, will excite a K-H instability, which will eventu- 
ally lead to the development of turbulence. A summary 
of the parameters we use is presented in Table ^ 

4. RESULTS 

4.1. Evolution of the mixing layer 

Our models, as described in the previous section, are 
dynamical and include a continuous transition between 
the hot and the warm media. In order to compare the re- 
sult of our calculations with previous models of turbulent 
mixing layers (SSB93) we have to define what a region 
of "intermediate temperature" is, in our computational 
box. We adopted a threshold of 20% departure from any 
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TABLE 1 

Parameters (initial conditions) of the runs used". 



iviouei 


l °g(J-hot) 


n-hot 


v hot 


D 

J-> 2 ,warm 


Grid size 




[KJ 


[cm -3 ] 


[km s ] 


[/XLrJ 




144-Th6-V50-B0 


6 


i x icr 3 


50 





144 3 


144-Th6-V200-B0 


6 


i x icr 3 


200 





144 3 


256-Th6-V200-B0 


6 


1 x 10~ 3 


200 





256 3 


LX-Th6-V200-B0 


6 


l x icr 3 


200 





256 x 144 2 


144-Th7-V50-B0 


7 


l x icr 4 


50 





144 3 


144-Th7-V200-B0 


7 


l x icr 4 


20 





144 3 


144-Th6-V50-Bl 


6 


2 x icr 3 


50 


~ 2 


144 3 


144-Th6-V200-Bl 


6 


2 x icr 3 


200 


~ 2 


144 3 


144-Th7-V50-Bl 


7 


2 x 10~ 4 


50 


~ 2 


144 3 


144-Th7-V200-Bl 


7 


2 x icr 4 


200 


~ 2 


144 3 



a All the runs start with a temperature on the warm side of T warm = 10 4 K and a hydrogen density of n warm = 0.1 cm 3 . The box size 
in the vertical direction corresponds to 10 pc. 

b All models were run either without cooling or with collisional ionization equilibrium cooling (see text for details). An additional prefix 
"NC" (no cooling) or "EC" (equilibrium cooling) to the run name is added respectively. 



of the nominal temperatures for the hot and warm media. 
That is, we consider in the transition region all material 
above 1.2 x 10 4 K, and below 8 x 10 5 K for the runs 
where T hot = 1 x 10 6 K; or bellow 8 x 10 6 K for the runs 
with Thot = 1 X 10 7 K. We evolve the initial conditions 
for ~ 20,000 time-steps for all the cases with 144 3 reso- 
lution. In Figures |5J and [21 we show the time evolution 
of all the purely hydrodynamical, and magnetized runs, 
respectively. 

We find in general, that the volume average temper- 
ature within the mixing layer depends strongly on the 
threshold used to define the intermediate temperature 
region. The density weighted temperature is more robust 
against the particular choice of threshold. And, since the 
emission observable is also density weighted, it is a better 
measure of the properties of mixing layers, and we will 
use it to interpret our results. Keeping in mind that the 
temperature in BF90 and SSB93 is a simple volume aver- 
age, and that denser regions correspond to lower temper- 
atures, our density weighted temperature will in general 
be lower. In our approach we do not assume a particular 
value for the efficiencies "77's, and/or £" , instead they 
are being implicitly calculated. In principle they can be 
measured from the simulations, however, it would require 
the choice of an arbitrary temperature threshold to ob- 
tain a mean (volume average) temperature, which turns 
out to be very sensitive to such threshold. Thus, it is 
impractical to compare the models through estimates of 
the efficiencies 

4.1.1. Formation of the mixing layer: early stages 

Using a density weighted average, our calculations 
show a relatively cold boundary layer (of a few 10 4 K) 
at the early stages of formation. Due to the fact that 
the mass is initially on the warm side, this indicates a 
somewhat inefficient mixing at such early times. 

Figures and show how the mixing layer devel- 
ops. When we do not include cooling, the thickness of 
the boundary layer increases monotonically. The growth 
rate is faster when the shear velocity is larger, as ex- 
pected for the K-H instability. For the cases in which 
cooling is included, we see a very different evolution, with 



Average Temperature 



1 o 7 r 



NC-**M-ThS-V - 
NC-144-Th6-V/0C-R0.0 ■ 
NC-144-Th7-V 50-RO.O - 



EC-144-Th7-V 




1 o 4 

0.0 0.2 0.4 0.6 0.8 1.0 1.2 
Time [Myr] 

H column density 




0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 
Time [Myr] 

Fig. 2. — Time evolution of the mixing layers, for the purely 
hydrodynamical runs. At the top (panel a) we present the den- 
sity weighted average (thick lines), and volume average (thin lines) 
temperature of the mixing layers (material with a departure of 20% 
from the nominal hot and warm temperatures). For reference we 
show the harmonic mean of the hot and warm temperatures (which 
corresponds to t in BF90, SSB93 with 5 = 1 (dotted) horizontal 
lines, see eq. 0). At the bottom (panel b) the mean hydrogen col- 
umn density in the mixing layer (both panels are averages over 
different lines of sight normal to the layer -along the x axis—). 
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Time [Myr] 




Time [Myr] 

Fig. 3. — Same as Fig. [3] for the cases with magnetic field 
(/3 ~ 1) on the warm side. 



the boundary layer shrinking with time at the beginning, 
and after some time starting to grow again. The small- 
est perturbations in the K-H instability have the fastest 
growth rate, but do not have the required energy to pull 
(enough) material into the mixing interface. Since the 
cooling is very effective, the transition layer is rather 
sharp. At this point, we mostly have warm gas being 
condensed at the interface, but mixing is not very effec- 
tive yet. This is demonstrated in Figure where we 
show cuts perpendicular to the mixing layer (XZ plane) 
in two of our simulations after t ~ 0.7 Myr. The cases 
with magnetic field do not show significant differences in 
the formation of the mixing layer when compared with 
the unmagnetized runs at early times (i < 0.7 Myr). 

In order to test the sensitivity of our simulations to 
resolution we ran a couple of cases in a 256 3 grid. The 
results are shown in Figure 

We found marginal dependence on numerical resolu- 
tion for the formation of the mixing layers. We will see 
however, that the mixing at later times is utterly domi- 
nated by the largest scales, when the modes that corre- 
spond to the longest wavelengths excited by the instabil- 
ity become apparent. 

4.1.2. Evolution of the mixing layers at later times 

Up until now, we have seen the formation of turbulent 
mixing layers by means of a K-H instability. However, 



we have not seen evidence of reaching a steady state, 
which is the original idea of the whole process. Since the 
largest modes in the K-H instability grow rather slowly, it 
is not practical to follow the evolution of all our models 
until steady state is achieved. Moreover, in a qualita- 
tive way, we have observed similar behavior for the runs 
that started with a larger temperature (T^ ot = 10 7 K ), 
but with the additional difficulty that we require finer 
time stepping for such cases, making it very difficult to 
cover a large span in time. To study the long-term evo- 
lution of the layer, we followed the the fastest evolving 
model. That is, the unmagnetized, T^ ot = 10 6 K, and 
v t = 200 km s _1 model. A related problem in our 144 3 
resolution runs was that by the time the mixing layer 
was reaching its asymptotic state (after ~ 1.5 Myr), it 
was getting very close to the boundary of the computa- 
tional box. To overcome this problem we extended the 
dimensions of the box to 256 cells in the x direction, and 
we placed the division between hot and warm media off- 
center. The new initial conditions ( "LX-Th6-V200-B0" 
in Table ^| have only 1 / 4th of the volume filled with 
warm gas, and the remainder with hot gas. This is be- 
cause we have much more heat capacity in the warm gas 
(due to the density contrast), and because our simula- 
tions did not have enough hot gas to provide a steady 
state. Figure El illustrates the evolution of this extended 
model, which shows how the mixing layer broadens after 
the first signs of large modes of the instability appear 
(see Figures 0Js and^Jf). 

We also show snapshots of cuts perpendicular to the 
layer, after 2.3 Myr in Figure [7| This figure highlights 
how a fully operational K-H instability effectively mixes 
the two media and provides a much larger interface zone, 
compared with the times when only small scale modes are 
present. Even after ~ 2.5 Myr, in the longer box (256 
cells in the x direction), the simulations have not reached 
a steady state. Nonetheless, our results provide impor- 
tant insights on how these mixing layers will develop and 
change over time. 

4.2. Observational Diagnostics for Mixing Layers 

Notably, the time evolution we observe for these mix- 
ing layers are associated with a significant evolution of 
spectral line diagnostics, most importantly the column 
densities and ratios of highly ionized species. These mea- 
surements depend not only on the metallicity, density 
contrast, and relative velocities of the mixing zones, but 
also on the evolutionary state of the layer. As a first step 
towards examining the evolution, we have calculated the 
column densities of C IV, N V, and O VI ions through our 
simulated layers as a function of time. Using the tem- 
perature and density profile along synthetic lines-of-sight 
(LOSs) through the data cube, chosen perpendicularly to 
the boundary layer (in the direction of the x axis). These 
column densities were computed under the assumption 
of solar metallicity, an d collisional ionization equilibrium 
ijBeniamin et alJl2001|) . Figures |H1 and [!|] show the range 
of values of the C IV/ O VI, and N V/O VI ratios along 
with other mo dels of high ion production in the literature 
(compiled bv llndebetouw fc Shulll2004albft . The central 
points are the mean C IV/O VI, and N V/O VI values 
averaged over all the possible LOSs, at a given point in 
time. The error bars were obtained at each time, with 
the dispersion of column densities for the various LOSs. 
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a) Time: 0.000 Myr Hydrogen Density [cm 

0J50 




b) Time: 0.000 Myr Temperature [K] 

I.00B+07P 1- 




c) Time: 0.707 Myr Hydrogen Density [cm 3 ] 

0.150 



d) Time: 0.707 Myr Temperature 




0.06b - 



0.02S - 




e) Time: 0.704 Myr Hydrogen Density [cm 3 ] 

0.150 



f) Time: 0.704 Myr Temperature 





Fig. 4. — Selected snapshots of cuts in the XZ plane from our simulations, showing density and temperature maps, along with the 
velocity field. The largest arrow in the velocity representation correspond to a magnitude of ~ 200 km s -1 The upper two figures (panels 
a, b) correspond show the initial conditions for vt = 200 km s — 1 , and T^ ot = 10 6 K. In the middle (panels c, d) we show the evolution 
after t ~ 0.7 Myr, without including cooling. At the bottom (panels e, fj we show the result after approximately the same time, but with 
the cooling included. We can see how the thickness of the mixing layer grows from the initial conditions where no cooling is present. At 
the same time, we can see the dramatic effect of cooling showing a rather sharp transition zone, and also how the largest wavelength modes 
of the K-H instability start to develop. 



The time domain was splitted evenly into 10 bins from 
t = to the longest time available for each run (the final 
times can be checked in Figures |3 El and 01 . 

The line ratios and column densities from our simula- 
tions should not be interpreted too literally, but rather 
as a guide to obtain insight into the evolution of mixing 
layers. In particular, the assumption of collisional ion 
fraction equilibrium limits the accuracy of line ratios. 
As pointed out by SSB93, mixing layers can be quite far 



from equilibrium. The overall dynamics and time evolu- 
tion of mixing layers is not likely to change dramatically 
by non-equilibrium effects. However, observational diag- 
nostics such as line ratios can be affected significantly. A 
study of non-equilibrium cooling models with more em- 
phasis on the observational implications is necessary and 
we plan to provide it in a following paper. 

It is interesting that a clear difference is present be- 
tween the runs that include cooling and those that do 
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Fig. 5. — Test of resolution convergence for a two set of simula- 
tions, one without cooling, and one with equilibrium cooling, both 
correspond to T^ ot = 10 6 K and vt = 200 KM s — 1 . In panel a) we 
show the density weighted average temperature (thick lines), and 
the volume average temperature (thin lines). Panel 6), shows the 
average H column density in the mixing layer. 



Fig. 6. — Further evolution of the runs corresponding to T^ ot = 
10 6 K, and v z = 200 m s — 1 . The top panel (panel a) with the 
average temperature, (thin lines for volume average, and thick lines 
for the density weighted.) At the bottom (panel b) we see that for 
the case with equilibrium cooling, after the initial stage, the the 
mixing layer starts to grow. 



not. However, note that except for the "LX-Th6-V200" 
runs, the mixing layer is still in a very early stage of for- 
mation. It is remarkable that the models without cooling 
show similarity to the results of SSB93. Certainly their 
models did not account for the dynamical effect of the 
cooling: their model provides the intermediate temper- 
atures, but the mixing is merely hydrodynamical. Con- 
densation or evaporating flows were not considered. At 
the same time, their model also assumed fully developed 
turbulence, and our set of runs have not achieved the 
corresponding stationary state. 

With new models, observations may be able to provide 
insights about the time-dependence of the mixing. An- 
other thing to notice is the larger scatter for the cases 
that have been run for longer times. This scatter can 
be explained noting that for longer timescales we have 
a more dynamical (turbulent) picture, with many struc- 
tures moving in and out of a particular LOS. Actually, 
the values of the C IV/O VI, N V/O VI ratios (and the 
rather large scatter), are somewhat consistent with th e 
observations presented bv llndebetouw fc Shulll (|2004b). 
Is is also evident that in the magnetized cases the K-H 
instability (and therefore the development of turbulence 
at the boundary) is delayed. This can be seen from the 



smaller scatter of points in Figure El compared to those 
in [HJ However, magnetic reconnection should not allow 
the magnetic field t o form knots and prevent t urbulent 
mixing motions (see lLazarian fc Vishniac|ll999l) . 

It remains difficult to reconcile the models with high 
ion col umn densities. It has been pointed out (e.g . 
SSB93: ISavage et al.l200.3l: llndebetouw fc Shull2004alH) . 
that to explain the typical column densities observed,the 
LOS must pass through several mixing layers (sometimes 
as many as a hundred). Indeed, several interfaces are 
likely to be blended for long lines-of-sight (as in the case 
of observations in the halo of the galaxy). However, in 
many cases fairly smooth and symmetrical line profiles, 
with little centroid dispersions are observed. This evi- 
dence indicates that merely summing over a large num- 
ber of interfaces may not explain the observed column 
densities. In Table [2] we present the column densities of 
an arbitrary artificial LOS after 0.7 Myr for all the mod- 
els with Th t = 10 6 K. For comparison we include the 
values that correspond to the initial conditions. 

In Figure 1101 we show the mean column densities of 
C IV, N V, and O VI (for synthetic LOSs normal to the 
interface), as a function of time for the more evolved 
runs. Average values of column densities in the halo 
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Fig. 7.— Snapshots of cuts in the XZ plane for "NC-L2X-Th6-V200-B0.0" after t ~ 2.3 Myr. Notice how the K-H instability has indeed 
produced turbulence at the boundary of the two media, increasing the thickness of the zone with intermediate temperatures (T ~ 10 5 K). 
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Ion Column Densities 


* AT EARLY TIMES. 






Model 




CIV 






N V 




O VI 




t = Myr. 


t = 0.7 Myr. 


t = 


= Myr. 


t = 0.7 Myr. 


( = Myr. 


( = 0.7 Myr. 


NC-144-Th6-V50-B0.0 
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12.06 




11.19 


11.20 


11.89 
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Fig. 8. — Comparison of line ratios for various mod- 
els in the literature, and our calculations (modified from 
llndebet ouw & Shull ( 2004a The triangles correspond to radia- 
tive cooling from a Galactic fountain model ( Shapiro & Bcniamin 
1993; Bcniamin & Shapiro 1993). The "x's" are conductive heating 
and e vaporation mod els of planar clouds jB ochringcr & Hartquist 
1981 IBorkowski. Balbus fc Fristroml [19901) Cooling of SNR 
shells are represented by "+'s" ISlavin fc Ccod 119931 ISheltonl 



19871 [borkowski, Bal bus fc Fristro: 



ds jBoehi 
119901). 



Cooling of SNR 



[1993). The ion ratios in c ollisional ionization equilibrium from 
Sutherland & Dopita (1993) arc in solid line (this corresponds to 
a range in temperatures, higher at the bottom center and lower at 
the upper right). SSB83 result for turbulent mixing layers are in 
asterisks. Our calculations are represented with circles, the color 
coding can be found in the legend of the plot. These last corre- 
spond to the average of LOSs normal to the boundary, for a given 
point in time each. The error bars were obtained considering the 
variability with the different LOSs. 



Fig. 9. — Same as Fig. |H1 with our simulations that include 
magnetic field (/3 ~ 1). 



High ions 




of the Galaxy observed with FU SE are ((Savage et alJ 
12003 llndebetouw fc Shull I2004H1 : logiV[C IV] ~ 14.3, 
logiV[N V] ~ 13.7, and log7V[0 VI] - 14.5. A survey 
with th e Copern icus satellite in the Galactic plane by 
Uenkinsl (|1978albf) . suggested that individual O VI com- 
ponents have column densities closer to log AT [O VI] ~ 
13, which is similar to observati ons with FUSE of very 
nearby stars ijOeeerle et alJl200"5j) . The column densities 
obtained in the model that ran for the longest time (at 
t ~ 2.5 Myr), are slightly larger than those predicted 
by SSB93, but would still require too many interfaces 
to explain observations. However, as we have discused 
earlier, the simulations were stopped before fully reach- 
ing a stationary state, and the thickness of intermediate 
temperature zone was still increasing with time (see for 
instance Figure m>). The results presented here suggest 
that not only the number of mixing layers, but the time 
available for turbulence to develop, are factors to con- 
sider for the proper interpretation of observations. 

5. SUMMARY 

Turbulent mixing layers at the interfaces of the hot 
(T - 10 6 ~ 7 K) and warm (T - 10 4 K) media in the ISM 



0.0 0.5 1.0 1.5 7.0 2.5 

Time [Myr] 

Fig. 10. — Time evolution of the column density of different high 
ions for a line of sight perpendicular to the interface. The initial 
conditions correspond to T hot = 10 e K and vt = 200 km s — ^ .Thin 
lines correspond to the case without cooling, thick lines to equilib- 
rium cooling. 



can be produced by high shear, via a Kelvin-Hclmholtz 
instability. We use a MHD code with radiative cooling, 
to model the formation and evolution of turbulent mixing 
layers through this instability. 

The introduction of cooling in to the dynamical evolu- 
tion was found to produce dramatic differences compared 
with the same models ran without cooling. For instance, 
the deposition of momenta from hot gas condensing onto 
the mixing layer modifies the growth-rate of the K-H in- 
stability, making it develop faster than the same case run 
without cooling. 

The low density of the hot medium caused the growth 
rate of the large scale modes of the instability to be 
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slow. At early stages of the formation of the mixing 
layer (i < 1 Myr) the rapid cooling of material at in- 
termediate temperatures (T ~ 10 5 K) makes for a sharp 
transition between the hot and warm media. 

At later times, the instability excites large scale fluctu- 
ations that are powerful enough to provide more efficient 
mixing, and the transition zone broadens. We see evi- 
dence that for typical ISM conditions (7\ot ~ 10 6 K and 
T W arm ~ 10 4 K, with a shear velocity of 200 km s _1 ), 
the mixing layer is still growing after 2.5Myr. 

We included a dynamically important magnetic field 
(~ 2 fiG) in the warm phase for some runs, and its effect 
was found to be minimal at earlier times < 0.7 Myr. At 
later times the magnetic field inhibits/delays the devel- 
opment of the K-H instability, and thus the turbulent 
mixing. 

Assuming solar metallicities, and collisional ionization 
equilibrium fractions, we estimated column densities and 
line ratios of highly ionized species (C IV, N V, and 
O VI) from our simulations. We compared our results 
with previous models (SSB93) of turbulent mixing layers 
and found similar column densities. The C IV/O VI, and 



N V/O VI mean line ratios are slightly different, with a 
lower C IV/O VI. We also found a significant scatter be- 
tween different lincs-of-sight, in particular for the models 
evolved to longer times. 

We were not able to fully reach a stationary state in 
our simulations, however our results suggest that given 
more time (> 2.5 Myr for typical conditions), the col- 
umn density of high ions should be significantly larger 
than previous models, partially alleviating the problem 
of number of layers required to explain observations. 
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